Multiparticle ring exchange in the Wigner glass and its possible relevance to strongly 
interacting two-dimensional electron systems in the presence of disorder 
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We consider a two-dimensional electron or hole system at zero temperature and low carrier densities, 
where the long-range Coulomb interactions dominate over the kinetic energy. In this limit the clean 
system will form a Wigner crystal. Non-trivial quantum mechanical corrections to the classical 
ground state lead to multiparticle exchange processes that can be expressed as an effective spin 
Hamiltonian involving competing interactions. Disorder will destroy the Wigner crystal on large 
length scales, and the resulting state is called a Wigner glass. The notion of multiparticle exchange 
processes is still applicable in the Wigner glass, but the exchange frequencies now follow a random 
distribution. We compute the exchange frequencies for a large number of relevant exchange processes 
in the Wigner crystal, and the frequency distributions for some important processes in the Wigner 
glass. The resulting effective low energy spin Hamiltonian should be the starting point of an analysis 
of the possible ground state phases and quantum phase transitions between them. We find that 
disorder plays a crucial role and speculate on a possible zero temperature phase diagram. 
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I. INTRODUCTION 

In recent years two-dimensional electron or hole sys=. 
terns with very low densities are intensively studied.EJ 
Such systems can be generated at the interface of gal- 
lium arsenide heterostructures or silicon metal-oxide- 
semiconductor field effect transistors, and more recently 
also in organic Cgo and polyacene films .el These materials 
provide an excellent environment to study the effects of 
strong electron-electron interactions and disorder. Qiic 
example is the unexpected metal-insulator transition.tH 

We consider two-dimensional electron or hole systems 
at zero temperature and zero magnetic field. In the ab- 
sence of disorder, it is known that the system will form 
a Wigner crystal in the limit of very low densities, where 
the non-trivial correlations can be-,described in terms of 
multiparticle exchange processes .Era The exchange fre- 
quencies then determine the magnetic Hamiltonian. A 
calculation of the exchange frequencies of a pure two- 
dimensional Wigner crystal was pioneered by Roger U 

Although conceptually important, the pure Wigner 
crystal cannot be realized in the systems mentioned 
above, due to disorder n A measure of disorder is the 
Drude conductance at an intermediate temperature scale 
at which the resistivity is relatively flat as a function 
of temperature, and the dominant contribution is from 
impurity scattering. At low densities, the measured 
Drude conductances are of order e 2 /h, indicating the 
importance of disorder. We consider this intermediate- 
temperature conductance as a tuning parameter for the 
quantum phase transitions to be discussed, not the 
asymptotic low temperature conductance. This charac- 



terization of the tuning parameter is important because, 
even for a pure system, the conductance at a 2D quantum 
critical point can be of order e 2 /hB 

It is also known that even an arbitrarily small amount 
of disorder will destroy the long-range order of the 
Wigner lattice.Q On short length scales, however, the lat- 
tice will remain unaffected by weak disorder, so that the 
notion of the multiparticle exchange is still valid. Strong 
disorder will compromise the crystalline order even on 
length scales comparable to the lattice spacing. Nonethe- 
less, the multiparticle exchange picture depends only on 
the existence of a rigid ground state in the classical limit 
(that is in the low density limit), which can be assumed 
to hold for any disorder strength. The exchange frequen- 
cies will, of course, follow a random distribution in the 
presence of disorder. _ 

In a previous paperO on the metal-insulator transition, 
we calculated a set of relevant exchange frequencies for 
the clean Wigner crystal within the many dimensional 
WKB approximation.!!^ This allowed us to conjecture a 
possible phase diagram in the ground state. The purpose 
of the present paper is to extend this calculation to the 
random distribution of exchange frequencies, necessarily 
caused by disorder in realistic situations. The resulting 
random and competing magnetic Hamiltonian should be 
an important ingredient in determining the phase dia- 
gram of this correlated complex system. A recent numer- 
ical calculation of exchange constants in a clean Wigner 
crystal is also available .Eil 



A. Wigner crystal and Wigner glass 

A two-dimensional electron system with carrier density 
n s is characterized by the dimcnsionless parameter 



- 1 =a B { m i s ) 1 / 2 , 



(1) 



which is a measure of quantum fluctuations; larger r s im- 
plies smaller quantum fluctuations. Here as — h e/m*e 2 
is the effective Bohr radius; m* is the effective mass, and 
e is the background dielectric constant. Thus, r s is the 
mean spacing between the carriers, in units of the Bohr 
radius. In a dilute system, where r s is large, we expect 
the ground state to be determined by the electrostatic re- 
pulsion between the electrons. In the absence of disorder, 
the classical ground state that minimizes the potential 
energy is a triangular lattice, the Wigner crystal. 

The crystalline state can be approximately described 
in terms of single-particle wavefunctions that locally re- 
semble harmonic oscillator wavefunctions. The spatial 
extent of these wavefunctions, Ar, depends on the oscil- 



lator frequency as Ar 



-1/2 



where w? is determined 



by the second spatial derivative of the electrostatic po- 
tential. A dimensional analysis yields u)$ ~ r s , so 
that Ar/r s ~ r s , and the system becomes increas- 
ingly classical as r s — > oo (n s — ► 0). At low densities, we 
can therefore systematically expand around the classical 
limit. 

As the density increases, or r s decreases, the Wigner 
crystal will melt at zero temperature. The melting tran- 
sition in (d + l)-dimension, where d > 1, is likely to be 
discontinuous from Landau theory formulated in terms 
of the ground state energy, which must be a unique funp. 
tional, E[p(r)], of the density, p(r), of the electron gas.Ej 
For a crystalline state, we can write 



(P( r )> = Pa+Yl P Ge ' 

G^O 



G : 



(2) 



where po is the average density and G's are the reciprocal 
lattice vectors of the crystal. In mean field theory, we can 
consider the ground state energy to be simply a function 
of the order parameters pc- Thus, the energy can be 
expanded as 



E = E[ Po ]+ l -Y, a G 



PG\ 
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The quadratic term is chosen to be 

a G = a(r^-r s ) + a / (G 2 -k 2 ) 2 , 



(3) 



(4) 



where a and a' are positive constants and ko fixes the 
length of the reciprocal lattice vectors of the crystal. For 
simplicity, u 3 and 114 were chosen to be momentum in- 
dependent, but functions of r s . On a triangular lattice, 
the cubic term is allowed by symmetry, hence the transi- 
tion to the crystalline state is discontinuous in the order 
parameter, or "first order" . 

Consider now the regime of the phase diagram for 
r s > r c s and weak disorder. We can prove that no matter 
how weak the disorder is the crystal falls apart at the 
macroscopic scale. It is sufficient to consider the limit 
r s 3> 1, because quantum fluctuations can only destabi- 
lize the crystal further. We can now apply the famous 
Imry-Ma-LarkinB argument. The gain in the pinning en- 
ergy due to disorder is proportional to L d / 2 , whereas the 
cost in the elastic energy of the crystal is L d ~ 2 , where L is 
the linear dimension of the sample and d is the space di- 
mensionality. Thus, for d < 4, the pinning energy wins, 
and the crystal is destroyed for arbitrarily small disor- 
der. Even if the crystal is disordered in the conventional 
sense, it still leaves open the possibility of a power-law 
ordered stateJia but this is now proven not to be possi- 
ble in d — 2.E-] The density-density correlation function 
fallsioff exponentially with a correlation length, £rj , given 

b y y 



£d > R a exp[cy/]n(Ra/a)]. 



(5) 



where R a is the length at which the displacement of the 
lattice becomes of the order of the lattice spacing a. Pre- 
cise calculations of the positive constant c, i? a , or the 
prefactor are not known for the Wigner crystal. Nonethe- 
less, £d is likely to be a large length in the limit of weak 
disorder, and it is safe to assume that short-range crys- 
talline correlations will survive. 

In d = 2, the lack of crystalline order, or even a power- 
law crystalline order, in the presence of disorder, does 
not allow us to argue for a distinct state of matter distin- 
guished by its special features with respect to the transla- 
tional degrees of freedom. From this perspective, one can 
continuously connect the liquid state and the amorphous 
crystalline state by moving into the disorder plane. Thus, 
in d = 2, the global symmetries that can be truly broken 
in the presence of disorder are the spin rotational invari- 
ance, 5, the time reversal invariance T, and the gauge 
invariance f (1). These symmetries can still label many 
distinct states of matter. For a related perspective on the 
problem of a pinned Wigner crystal in a magnetic field, 
see Ref. 16. We note that in d — 3 a power-law ordered 



Wigner glass can exist as a distinct state of matter. 



B. Magnetism: pure system 

In discussing the magnetism of the insulating Wigner 
crystal, we shall ignore anharmonicities of the zero point 
phonon degrees of freedom, which may merely renor- 
malize the exchange constants. The low lying magnetic 



Hamiltonian is due to tunneling of electrons between 
the lattice sites and can be expressed in terms of the 
p-particle cyclic permutaion operators P\... p . Thus, 

H = J *Y. ( Pi - 2 + p_1 ) - J3 X! ( Pi > 2 - 3 + p_1 ) 
~ A 

+ j 4 Y. (pi-a + p- 1 )-^ J2 ( p i.5+^ 1 ) 

+ J 6 Y. ( p i-s + P~ 1 ) + ---- (6) 
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FIG. 1. Schematic setup of the system in which a two- 
dimensional electron or hole gas is generated. 



The sums are over the permutations shown in this equa- 
tion. There is a theorem due to Herring and Thouless 
that exchanges involving even number of fermions are an- 
tiferromagnetic, and those involving odd number of par- 
ticles are ferromagnetic. a We shall follow the convention 
that the J's are all positive. 

A tractable method for calculating the exchange con- 
stants J2, J3, ... is the instanton (or the many dimen- 
sional WKB) method. It will be shown that J p is 



J P = A p h^ (A.) 



1/2 



-S p /h 



(J) 



where S p is the value of the Euclidean action along 
the minimal action path that exchanges p electrons. 
The quantity ujq is the characteristic attempt frequency, 
which can be estimated from the phonon spectrum of 
the lattice. The prefactor A p is of order unity, and the 
Eq. (0) holds as long as % > 1. 

The cyclic permutation operators can be expressed in 
terms of the spin operators using the Dirac identity P12 = 
-^ + 2Si • S2 and the spin Hamiltonian is 



H 



Jnn / &i 



Jn 



z^ 



Oi ' Ji 



(8) 



The first term in Eq. (0) is sum over distinct nearest 
neighbors, the second is over distinct next nearest neigh- 
bors, and so on. Here, J nn — 4J 2 + 5 J4 — 4 J 3 + • • ■ and 
Jnnn = J a + ■ • '■ In general, this is a highly competing 
magnetic Hamiltonian. On a regular triangular lattice a 
model containing exchanges upto J5 has been studied by 
various approximate analytical and numerical fait^ size 
(maximum of 36 sites) diagonalization methods .Ota The 
picture that has emerged is rather complex containing a 
number of broken symmetry states: a ferromagnetic, a 
three sublattice Neel, a four sublattice Neel, and a long 
wavelength spiral states. In addition, on the basis of nu- 
merical work, it has been argued that a sizeable region 
of the phase diagram consists of a spin liquid state, with 
short ranged correlations, spin gap, and no broken trans- 
lational and spin rotational symmetries. 



C. Magnetism: disordered system 

In the presence of disorder, the picture should change 
substantially. The system is no longer described by a 
regular triangular lattice and will instead distort into a 
random lattice, with the sites dominantly determined by 
the pinning defects. Those properties of the pure system 
that are specific to a triangular lattice will no longer hold. 
For example, none of the antiferromagnetic states, which 
depend delicately on the regular lattice structure can be 
the true ground states. More fundamentally, there is 
no longer an argument that 3-particle exchange is larger 
than the 2-particle exchange, rather the opposite could 
hold, as we shall see. To explore the effect of disorder, 
we calculate the multiparticle exchange processes in a 
disordered system whose low energy magnetic Hamilto- 
nian can be formulated as in the pure system but with a 
random distribution of exchange constants. 

Leaving aside the gauge symmetry, the symmetries 
that are allowed to be broken in a disordered system are 
the spin rotational invariance (S) and the time reversal 
invariance (T). The phases that are potentially impor- 
tant are a T-broken metal, a T-broken insulator, a S and 
T broken spin-glass, a disordered ferromagnet, and a dis- 
ordered antiferromagnet. Since a disordered system does 
not respect translational invariance, no further subclas- 
sification according to broken translational symmetry is 
possible. It is clear, however, that the regime close to the 
crystalline phase of the pure system will be marked by 
strong short-ranged crystalline order. Generically, dis- 
order necessarily renders all quantum phase transitions 
between these states continuous, and thus the phase di- 
agram is rife with quantum critical points and lines. 



D. The Model 

The systems of experimental interest differ consider- 
ably, but they can be schematized as shown in Fig. |lj. 
The carriers themselves are confined to an inversion layer 
or a quantum well with a width of the order of ~ IOOA. 
A buffer of several hundred A separates the carrier plane 



from the doping layer, which contains impurities in the 
form of oppositely charged ions that provide the carriers. 
We will use the language appropriate to the electron- 
doped case for the sake of clarity; the hole-doped case 
can be treated identically. 

Let us denote the coordinates of the TV electrons by rj, 
and those of the Ni mp positively charged impurities by 
imp ^ ^y e -vvill treat the carriers as being exactly confined 
to the xy plane, so that r^ = (#,-, j/j, 0), which means that 
we neglect the finite spread of the wavefunction in the di- 
rection perpendicular to the plane. This spread leads to a 
softening of the Coulomb potential at distances compara- 
ble with or smaller than the effective Bohr radius as ■ In 
the dilute limit consider here, r s 3> 1, the many-particle 
wavefunction will be negligibly small in those regions of 
coordinate space where two or more electrons come close 
enough to each other to "feel" this softer potential. 

We assume the only source of disorder is provided by 
the impurity ions in the doping layer, which is separated 
by a distance d from the carriers. We will consider the 
following model for the impurity distribution: the thick- 
ness of the impurity layer is taken to be zero, so that the 
impurities are exactly confined to the plane z = d. We 
have also considered a second model in which we assumed 
the impurity layer to have a finite thickness, taken to be 
equal to the separation d from the electron gas. Since the 
results are very similar, we shall not report them here. 

Within the doping layer the impurities are randomly 
distributed. The Hamiltonian is 



N [ 2 i-1 N imp 

»=1 l 3=1 3=1 



imp , 



Z = E (-l) P E/rfRG(R )£ ;PR,P £ ;/3) (12) 

P&S N a 

where the first sum is over all N\ permutations PR of 
the electron coordinates, and (— l) p is the sign of the 
permutation. The propagator is defined as 



G(R.i,a 1 ;R 2 ,£ 2 ',T) 



(Ri,*i\e-™ \R2,a 



(13) 



Here 5 a a is a product of N Kronecker delta symbols. 
Note that this definition of the propagator treats the 
electrons as distinguishable Boltzmann particles. Fermi 
statistics have been taken into account in the sum over 
boundary conditions in the partition function (|l2|). 



A. The Semiclassical Approximation 

The instanton method that we-.shall follow has been 
elegantly discussed by Coleman. Ej The imaginary time 
path integral for the propagator is (T T here denotes imag- 
inary time) 



G(Ri , g_ x ; R2 , g_ 2 ', T r ) = S ZiZi 



where the Euclidean action is 
^n /" Tt , [m* fdR 

s t R ]= L * -r * 



R(T T )=R 2 



R(0)=Ri 



DRe^ s[R1 



V(R) - V 



(14) 



(15) 



where 



. e 2 1 

v(r) = 7w 



(10) 



is the effective Coulomb potential, e being the dielectic 
constant of the environment and m* the effective mass 
of the carriers. 



II. THE MULTIPARTICLE EXCHANGE 
PICTURE 

It is useful to define the collective spatial and spin co- 
ordinates 

R= (n,r 2 ,---,rjv), g_= (cri,o- 2 ,---,o-jv). (11) 

Formally we can view R as the coordinate of a single par- 
ticle moving in a 2iV-dimensional space in the potential 
V(R); see Eq. | 

For Fermions, the partition function of the system is 
then 



The equilibrium potential energy Vq = minR,y(R) has 
been subtracted out for later convenience. The stationary 
path satisfies 



6S[K C ] ,d 2 TL c -,_ . n 

and the action for this path is 
Rf 



S[Rc] = / dR y/2m* (V(R) - V ). 



(16) 



(17) 



R, 



The Planck constant enters the action in this form only 
through the parameter r s ~ l/h . Therefore, the semi- 
classical caluclations described here are accurate in the 
low density limit, r s — * 00. 

The Gaussian quantum fluctuations around the sta- 
tionary path are taken into account by defining the fluc- 
tuation coordinates u(r) = R(r) — R c (t), in terms of 
which we expand the action to second order: 



f]JR e -^ R ] = F[R c 



e R 



SfRc 



(18) 



where 



F[R c ]=/Duexpj-i^ T Vu 



(r) • A(r)u(r) 



[detA 



-1/2 



(19) 



and we have assumed that the stationary path is unique. 
In cases where more than one stationary path exists, their 
contributions have to be summed. The differential oper- 
ator A is given for each path by 



A^ v {t) = -5fj, v m* 



d 2 V{K) 



dr 2 dRudR,, 



(20) 



R=R c (r) 



The determinant is defined in terms of the eigenvalues X v 
of A, subject to the boundary conditions u(0) = u(T) = 
0, as Act A = Y\ V K- 



B. Exchange Processes and the Instanton 
Approximation 

We will assume that there exists a definite configura- 
tion R of the N electrons that minimizes the electrostatic 
potential V(R). It is clear that this classical minimum 
is ATI-fold degenerate, since the potential energy is in- 
variant under any permutation of the electron coordi- 
nates. In the semiclassical limit, configurations where R 
is in the vicinity of one of these minima will contribute 
dominantly to the partition function. We will therefore 
construct stationary paths that begin and end at one of 
those minima. In particular, we define the instanton path 
Rinst(T) between the two minima at Ri and R2 such that 



Rinst(-oo) = Ri and R inst (+oo) 



R 2 



(21) 



and the equation of motion ( |16[ ) is satisfied. In the sim- 
plest case, that of Ri = R2, the instanton path is given 
by Ri nst (r) = Ri. In general Ri and R2 differ by a 
permutation of the electron coordinates r^, so that the 
instanton path describes a multiparticle exchange pro- 
cess. 

In the vicinity of a minimum, R, 



R(r) - R . 



e ±^( R )u , 



(22) 



where u is some constant vector, and fl is defined as the 
square root of the Hessian matrix, evaluated at R: 



_. 1 d 2 V(TL) 



njUR) = 



dR^dR^ 



(23) 



R=R 



Hence any deviations from the classical equilibrium con- 
figuration are localized on the imaginary time axis on a 
scale St ~ -j-, where uu 2 is some eigenvalue (not neces- 
sarily the smallest) of f2(R). In this sense, the instanton 
path will be localized around the location of the instan- 
ton, Tjnst) in imaginary time. On a coarse-grained time 
scale the exchange processes will therefore appear as in- 
stantaneous, independend events. 



The instanton path will be unique in most cases. An 
exception is the two-particle exchange, where the elec- 
trons can take two equivalent paths corresponding to 
clockwise and counterclockwise exchange. This merely 
results in a factor of two for the exchange frequency. 

The instanton formalism rests on the assumption, oc- 
casionally referred to as the dilute gas approximation, 
that the average distance At on the imaginary time axis 
between exchange processes within the same region of 
space exceeds the instanton duration St by several or- 
ders of magnitude. If we consider the propagator on a 
time scale T T that satisfies 

5t<^T t <^ At (24) 

we can make the following two crucial assumptions: 

1. Each time slice contains at most one instanton 
event. Processes with two or more instanton events 
in a single time slice are of second order in T T j At 
and therefore negligible. 

2. Instantons do not occur within a few instanton 
lengths of a time slice boundary. Again, pro- 
cesses that violate this assumption arc of order 
(T t /At)(St/T t ) = 5t/At and therefore negligible. 

Let us now evaluate the propagator within these approxi- 
mations. Since the Hamiltonian is independent of spin we 
will suppress the spin indices in our notation for the mo- 
ment. We also define a fluctuation coordinate u = R R, 
where R is by definition the particular minimum of V(R) 
that is closest to R. Thus we want to evaluate 

_ /.R(ti)=R 2 +u 2 

G(Ri+Ui;R 2 +u 2 ;T r )= / _ DRe^ s ' R ], 

JR(0)=Ri+ui 

(25) 

where the deviations Ui and U2 are by assumption 2 in 
the quadratic regime, so that we can expand the equation 
of motion to linear order in u. In other words, we are 
allowed to approximate 



o(r + u) ~n(R). 



(26) 



An approximate solution of the equation of motion dig ) 
that satisfies the boundary conditions 

R c (0) = Ri + Ui and R C (T T ) = R 2 + u 2 (27) 

is then 

R c (r) - e-^ /2 Ul + eS T - T ^< 2 u 2 + R inst (r - r ), 

(28) 

where tq is an arbitrary reference point between and t\ . 
The time derivative of each term is localized on a time 
scale St <C T t , and by assumption 2 above the overlap 
between the three terms is exponentially small. Hence 



the corrections arising from the nonlinearity of the equa- 
tion of motion are negligible. For the same reason the 
action associated with this path splits into three parts, 
which we write in an obvious notation as 



5[R C ] = S[ui] + S[u 2 ] + S 



inst- 



(29) 



With the results of the previous section the propagator 
is then 

G(Ri + ui;Ra + u 2 ;T r ) 

= F [Ri n8 t] exp (~ (5[uJ + S[u 2 ] + S inst )\ , (30) 

where we already incorporated the fact that with the 
approximation (U&j the prefactor (03) is independent of 
the Ui . Let us from now on write Ri, ls t = Rp and S- Uls t = 
Sp, where P G Sn labels the particular permutation that 
takes R2 into Ri, i.e. Ri = PR2. For later use we define 
the quantity Gp as the ratio between the propagator for 
a given instanton path and the corresponding propagator 
for the trivial path Ri nst (r) = R, which has S^st = 0. 
That is, 



G P := 



G(R + ui;PR + Pu 2 ;T T ) F[R P ] _ l3p 
= = = — ^=^- p. 1 r . 



G(R + ui;R + u 2 ;T T ) 



F[R] 



(31) 



Within the instanton approximation Gp is indepen- 
dent of the fluctuation coordinates. Due to permutation 
symmetry all terms in Gp are also independent of the 
particular choice of the minimum R. 



C. The Prefactor 

To evaluate the prefactor ( |l9|) we would have to find a 
complete set of eigenfunctions u n (r) that satisfy 



~m*5 tJllJ d 2 + V[j,v(t)] u nl ,(r) = A„u, iai (t) 



(32) 



with the boundary conditions u n (0) = u„(T T ) = 0. 
Here we used the shorthand notations d 2 — d 2 /dr 2 and 
V^(t) = d 2 V(R P (T))/dR,idRv. If we expand 



u(r) = ^c„u„(r) 



n=0 



the path integral over u is transformed into 

dc„. 



Du 



n 



(33) 



(34) 



We still have to account for the possibility that one 
of the eigenvalues of ( P0|) is less than or equal to zero. 
While it is easy to show by direct calculation that this 
is not the case for F [R] , a zero eigenvalue indeed exists 



for F [Rp] . As an eigenfunction we consider the time 
derivative of the instanton path itself: 



(35) 



u (t) := 


a ° _1 | Rp(r - 


To), 


where the normalization constant is given by 


a 2 ) = dr 
Jo 


d _ , s 


\Sp 

TO* 



■±- (36) 



The last identity follows from the equation of motion 
(|l6|), which can be integrated to give 



m 



dn c 
~dT 



V(R C ) - V , 



(37) 



which is just the Euclidean version of energy conserva- 
tion. It is straightforward to verify that Uo(r) satisfies 
the eigenvalue equation ( p2[ ) with eigenvalue Ao = 0. 
The boundary conditions u(0) = u(T T ) = are satis- 
fied within our approximations since Uo is exponentially 
localized. This function just describes the change in R c 
due to a change in the instanton position to, which is 
arbitrary within the limits < tq < T T . Hence a shift 
in the instanton position corresponds to a zero mode. 
This is the Goldstone mode associated with broken time 
translation symmetry in the presence of an instanton. 
The change in the path R(t) = Rp(r) + u(r) due to a 
change in the expansion coefficient cq can be related to 
a change in tq as follows: 



dR(T) 
dc 



Uq(t), 



byEq. @, while 

dER(r) _ dR P {T-T ) 



dm 



dr 



= -a u (r) 



(38) 



(39) 



by the definition of Uo, so that we have to replace the 
integration over cq by 




= T T 



Si 



2n~hm* 



(40) 



Ao is the lowest eigenvalue, since the corresponding 
eigenfunction is free of nodes. Hence all other eigenvalues 
must be positive. We now have 



FRl ' :7 W^J^(det'[-m^+ VAW (T)]) 



-1/2 



(41) 



where the prime indicates that the zero eigenvalue has 
to be omitted in the determinant. To summarize, the 
prefactor is given by 

1/2 



F [Rp] _ / S P (det[-m*d 2 + V^(0)] 

F[R] ='- 



2Trhm* \ det' [-m* d 2 + Y^ (r) 



(42) 



Let us assume that we scale all eigenmodes of the po- 
tential by the same factor g and simultaneously rescale 
the imaginary time variable by a factor g~ l . The factors 
of g in the determinants cancel in the numerator and in 
the denominator for each eigenvalue separately, and we 
know that the ratio of determinants cannot depend on the 
length T T of the time interval, except for exponentially 
small corrections. Hence the prefactor depends linearly 
on a characteristic frequency scale of V(R), and the ra- 
tio of determinants depends only on the relative values 
of the eigenfrequencies. We summarize these findings by 
writing Gp, defined in Eq. (pl|), as 



Gp — T T Ap luq 



2-kTi 



■isp 



(43) 



where, as stated above, ojq is a characteristic frequency, 
and the dimensionless factor Ap depends on the relative 
values of the eigenfrequencies during the exchange pro- 
cess. It seems reasonable to assume that Ap is roughly 
of order one. Although the prefactor is not expected to 
cause any drastic changes in our results, it is still interest- 
ing to determine the change in characteristic frequency 
with disorder. It is conceivable, for example, that disor- 
der would bring about a reduction in the phonon spec- 
trum, and this mechanism could lead to a suppression of 
exchange processes. 



D. The Exchange Hamiltonian 

Our goal in this section is a Hamiltonian description of 
the system in terms of multiparticle exchange operators. 
In the previous subsections we calculated the imaginary- 
time propagator on an intermediate time scale T T de- 
fined by the relation ( p4|). To apply this result, we split 
the partition function (|l2|) into M imaginary time slices, 
where M satisfies (3 = MT T . This requires us to sum over 
M — 1 intermediate configurations, so that the partition 
function reads 



^ = E(-d p E---E 

p 









dRi 



r/R 



M 



-i 



-M 



x G(Ri,o; 1 ;R 2 ,a 2 ;T r ) 

■■■G(Hm,<Lm\P^P<Li\Tt). 

We want to make use of the quantities 

G(R + u i; PR + Pu 2 ;T T ) 
Gp = — — = ==- 



G(R + ui;R + u 2 ;T T ) 



(44) 



(45) 



defined in Sec. 



II B 



5, which only depend on the permu- 
tation P, and are independent of the fluctuation coor- 
dinates Ui and the particular choice of the minimum R. 
To this end, we write the integration variables Ri in the 
form 



R i = P i (R + u i ), 



(46) 



where R is some minimum of V(R), and the permutation 
Pi is chosen such as to minimize the distance |R-j — PjR|. 
The integrals then have to be replaced with 



/ dRi -»E / dut, 



(47) 



where the sum is over all permutations, so that P^R cov- 
ers all minima of V(R), and the integration over u^ is by 
construction restricted to the vicinity of u^ = 0. The par- 
tition function then reads (dropping spin in the notation 
for now) 






r/u 



M 



(48) 



p Pi p, 

x G(Pi(R + u 1 );P 2 (R + u 2 );T T 

■ • • g(p m (R + u M );PPi(R + ui);T T 



Let us now introduce the transfer matrix T, defined by 
the relation 



G\P i (R + u i ),s i )P 3 -{R + xi j ),a j ',T Tj 

= (i,2i\f\j,2j) G{R + Ui;R + u f ,T T ). (49) 

Comparing this definition to Eq. ( j45| ) we can easily de- 
duce 

(i,£i\ T \j:°-j) = 5 g_ lS _ 2 G Pzj = E Gp (*'2»l *" b'>£j) ' 

p 

(50) 

where P^ = P^ l Pj, the permutation operators P' are 
defined to act only on the index i as P' \i, a) = \pi, a), and 
we made use of the orthogonality relation (i,g^\j,aj/ = 
SijS a a ■ Thus the transfer matrix is 



t = j2 G p p '- 



(51) 



Inserting Eq. (49) into the partition function (48), the 



latter will factorize into a fluctuation part and a tunnel- 
ing part: 

^oE(-d p E-EE-E 

P h *m 5-i S-m 

x (ii,2: 1 |T|i2,22) ■■■{iM,g: M \TP\ii,a 1 } 

= Z„^(-l) p ^( i ,a|T M P| i ,a), 



where the permutation operator P acts on both i and g^ 
as P\i,g) — \Pi,Pg), and 



Z {) = I dm-- / du M G(R + ui;R + u 2 ;T T ) 

•G(R + u M ;R + u i; T r ) (52) 



is the partition function for a 2./V-dimensional harmonic 
oscillator. 

The Gp are proportional to the lenght of a time slice 
T T = 0/M, see Eq. (|43|), with the exception of the iden- 
tical permutation P = 1, for which G\ — 1. We therefore 
define the exchange energies Jp = ^-Gp, which allows 
us to write 



rpM _ 



£ 

M 



j2j p P' 

Pj,\ 



M 



= exp 






in the zero-temperature limit, in which M = 0/T T 
The partition function now reads 



(53) 



oo. 



P\ 



12 



2Si • S 2 






1 /1Z-2 

2\ a \ 2 



(58) 



leading to a Heisenberg term, as one can easily check 
by direct calculation of the matrix elements. Any per- 
mutation can be written as a combination of elementary 
transpositions, and hence as a product of spin operators. 
In general these products can be reduced using operator 
identities such aalJ 



per , pc _ per I per , per _ i 

-"123 ' -"321 — r V2 ' r 23 ' r 31 1 ' 

per i per per per , per f)cr per per 

^1234 T ^4321 — ^12^34 "+" ^14^23 r 13 r 2A 



~Pl 3 



per _ i 

-"24 l ) 



(59) 



P ia [ P'jtl J 

(54) 

This is the desired representation in terms of permu- 
tation operators. The exchange energies are given by 
Eq. (|H) as 



Jp=A P hu; aX l—e 



-*s r 



(55) 



If we expand the exponential in a 
the orthogonality condition (i,g_\j, g_) - 
that all permutation operators P' in this expansion 
have to combine with P to the identical permutation: 
P[P!, ■■■P^P |», <r) = |t. Pa), or P = (P^)" 1 ■ • • (P[)-i as 
far as their action on i is concerned. We can thus elim- 
inate the sum over P and absorb the spin permutations 
and the sign factor into the exponential. Then the sum 
over i is redundant due to permutation symmetry and 
the partition function becomes 

Z = N\Z ]T (a\ exp J -0 ^(-1) P+1 J P P° \ \a) , 

(56) 

where P a acts on the spin variables as P a \a) — \Pg). 
This is the partition function for a pure spin Hamiltonian 



H, = J^(-1) P+1 JpP° 



(57) 



E. Generalized Heisenberg Model 

The spin permutation operators appearing in Eq. (|57 
can be rewritten in terms of Pauli spin operators. For 
example, if we denote by P° 2 the permutation operator 
that interchanges o\ and C2, 



etc. Keeping only the dominant two-, three-, and four- 
particle exchange processes, the spin Hamiltonian be- 
comes 



NN 



NNN 



H = (2 J 2 - 4J 3 + 2 J 4 ) Y^ s i ■ s i + 2J i Yl Si • Sj 



<y> 



<y> 



-4J 4 2_^ (Gijki + Gujk — Gikji) , 

<ijkl> 



(60) 



where NN indicates a sum over nearest-neighbor pairs, 
NNN a sum over next-nearest neighbors, and o is a sum 



power series, over all rhombi. G^ki — (S, • Sj)(Sjt ■ Si), where the 
= Sij implies vertices of the rhombus are labeled clockwise by the four 



indices. This Hamiltonian has been discussed in Sees. 



and IC 



[B 



III. NUMERICAL TECHNIQUES 
A. Calculation of the Action 



The action (fil) 



Sp = _ dR^2m* (V(R P ) ~ V ) 



(61) 



depends only on a single length scale, which can be fac- 
tored out. We define a dimensionless coordinate X = 
iR, where a is the lattice constant of the ordered Wiener 
crystal. The unit cell of the triangular lattice is of area 
A = -^a 2 , so that the density is 



1 
A 



V3 



(62) 



which we use as a definition for a in the pre senc e of dis- 
order. The parameter r s . 
expressed in terms of a as 



defined in Sec. I A, can be 



3 1 ' 4 a a 

r s = —== — -0.525 — . 

V27T a B CLE 



(63) 



We define the dimensionless action Sp by j-Sp 



r s Sp. Then 



PX 



S P =v _ dXdV(X P )-V Q 



where 






N I i-X 



2 = 1 I J — 1 ' Jl J = \ JS^ Aj 



(64) 



(65) 



is the dimensionless potential, Vq is its minimum value, 
and 77 is a numerical factor: 



-*(*) 



1/4 



1.952. 



(66) 



The classical path that minimizes the action has to be 
found numerically. Therefore we discretize the integral 
in (64) using the trapezoidal rule, which leads to 



M-X 

S — o / y l x m - X 'i 



V(Xi) - F + \MXn-i) " ^c 



(67) 



The displacement of the participating electrons from 
their equilibrium positions creates dipole perturbations, 
which are screened out after a distance of a few lattice 
spacings, even in the absence of conventional screening. 
We can therefore restrict the number of moving particles 
to a relatively small value iV m obiic and hold all other par- 
ticle coordinates fixed at their equilibrium values. De- 
tails on the errors due to the finite values of M and 
-^mobile can be found in Appendix |A[ Since these er- 
rors are of opposite sign, we believe that the total error 
for the action is no larger than 0.3% in the clean sys- 
tem. In order to keep the distances |X<+i — X»| approx- 
imately constant during the minimization process, the 
allowed variations in Xj are restricted to those satisfy- 
ing <JXj • (PX — X) = 0, thereby reducing the number 
of independent variables per time slice by one. Since ini- 
tial (i = 0) and final (i = M) conditions are held fixed, 
the action is a function of (27V — 1)(M — 1) independent 
variables in its discretized form. For calculations on the 
clean system we took M = 16 and iV mo biie — 80, de- 
pending on the particular exchange under consideration. 
The minimization thus involves around 2400 variables. 
We used a variable metric (quasi-Newton) algorithm.E2l 
Due to the long range nature of the Coulomb potential 
the sum in the expression ( pq ) for the potential energy 
converges very slowly, and is in fact only conditionally 
convergent. We therefore use the Ewald summation tech- 
nique, in which the summation over the long-range part 
of the Coulomb potential is carried out in Fourier space. 
To improve the speed of the computation, we tabulated 



the Ewald summation formulas on a 50 x 50 grid and 
calculated in-between values using bicubic interpolation. 
We explicitly checked that the interpolation procedure 
does not generate any errors comparable to the stated 
accuracy of the results. In this way a single minimiza- 
tion could be carried out in less than 10 minutes CPU 
time on a 400MHz Pentium II processor. 



B. The Prefactor 

We now turn to the numerical evaluation of the pre- 
factor (19), which we write in the form 



F[R C ] = f Du 



-s\u\ 



(68) 



dr u(r) z + u(t) • #(t)u(t) , 



where 



H„»(t) = — 



1 d 2 V(R) 



m* dRndR v 



(69) 



R=R c (r) 



is the Hessian matrix of the potential for the configura- 
tion at time r. We split the imaginary time axis into M 
intervals, so that 



T T = s Q > si > • ■ • > sm-x > s M = 0, 



(70) 



and approximate H (r) by a constant matrix Hi within a 
given time interval s, > r > s-j+i: 



H^(t) ~ {Hi)^ 



1 d 2 V(R) 



dR^dR v 



(71) 



R=R, 



where R^ = aXj are the point s of t he discretized instan- 
ton path determined in Sec. Ill A. The corresponding 
times Si can be calculated by inverting the equation of 
motion: 



Si+l 



dR 
m* 



m 

1/2 



:(V(R)-V ) 



-1/2 



(72) 



R»+i — Rj + IRj — R,--i 



We can then write the prefactor in the form 
F[R C ] ~ /"duiGi(0,ui;T T -si) 
x / du 2 G2(ui,u 2 ;si - s 2 ) 



■ Gm (um-i, 0; sm-x), 



(73) 



where 



Gi{ui,Ui +1 ;s) = / Du cxp i -— / dr 



du' 



u(r) ■ Hi u(r) 



(74) 



is simply the propagator of a multidimensional harmonic 
oscillator and can easily be calculated. We define or- 
thonormal eigenvectors e\, and eigenvalues oJf v satisfying 



Hi el = <4u K 



H c„ — LU it/ c„ (75) 

(note that u>i V can be imaginary) , in terms of which 




\ 1/2 
\{B lv { S )) x 



~(u ul + u 2 v2 ) cosh u) u s - u u iu v2 



m LUi V 



(76) 



?jsinhw,-,,s 



(77) 



The prefactor is then 



F[R C ]= I Y[B iv ( 

\ iv 



-1/2 



Si - Si-l . 



(detM)" 1/2 , (78) 



where the matrix M is defined as 



M% = 5 id (4„ + A- 1 ) - (Sij+iB^ + Sij-iBJJ) , 



A]„, = 



^r 



m uj iX 



'"" ^ ^ A "tanha, iA (*-*_!)' 



■^*v — 2^ e \p. e Xv 



m to l \ 



smhcji\(si - Si_i) 



(79) 



Numerical evaluation of the determinant is now 
straightforward. The e igenvalue that corresponds to the 
zero mode of Sec. [I C has to be omitted from the result 



(this eigenvalue will not be exactly zero here, due to the 
finite number of time slices). To this task, we replace 
H(t) by H{t) — A in Eq. (pq), and numerically search 
for the smallest value of A that satisfies 1/F(A) = 0. We 
then divide the determinant by this value. The method 
outlined here has been tested on the problem of tunnel- 
ing in a quartic potential in one dimension, which can be 
treated analytically; details can be found in Appendix M. 



C. The Disordered System 

For the disordered system we have to sample over a 
large number of impurity distributions. After placing the 
impurities onto random locations in systems with 48 to 
280 particles and periodic boundary conditions, we first 
minimize the potential energy of the classical electron 
configuration. No tabulation of the Ewald summation 
formulas was used in this minimization, since the clas- 
sical equilibrium configuration is very sensitive to nu- 
merical errors. We cannot exclude the possibility that 



the minimization procedure gets trapped in a metastable 
configuration in the presence of strong disorder. On av- 
erage, however, the properties of such a metastable state 
should be sufficiently similar to those of the true equi- 
librium state that our results will not be affected. For 
strong disorder, when the triangular lattice structure is 
compromised even on short length scales, we are also 
faced with the problem of identifying proper sets of near- 
est neighbors to participate in the exchange. This task 
is solved by a Delaunay triangulation of the electrons' 
equilibrium coordinates. For the subsequent minimiza- 
tion of the discretized action only the N mo bn e = 32 • • • 34 
particles closest to those participating in the exchange 
were allowed to move, with the remaining particles held 
fixed at their equilibrium positions. The number of time 
slices was reduced to M = 8, so that we have to mini- 
mize over approximately 500 independent variables. The 
minimization converges significantly slower than in the 
clean system, since the dependence of the action on the 
independent variables is less smooth. In the presence 
of strong disorder each minimization takes several min- 
utes to carry out. Typically we generated around 250 
impurity configurations, for each of which 8 exchange 
processes were chosen at random between sets of near- 
est neighbors anywhere on the lattice. We thus arrive at 
about 2000 sample values per data point. 



IV. RESULTS 

A. The Clean System 

Here we present results for a large number of exchange 
processes in the absence of disorder, including all those 
that are relevant at low densities. The exchange paths 
are shown schematically in Fig. |2j, and the correspond- 
ing values of the dimensionless action S n are listed in 
Table Q. Roughly speaking, the action depends both on 
the number of particles involved, and on the smoothness 
of the exchange paths. Kinks in the path are penalized, 
since they lead to intermediate configurations with high 
potential energy. This is also the reason for the relatively 
high value of S 2 ■ For the smoothest exchange paths with 
n > 8 the action increases roughly linear with n (see 
Fig. g). We have, approximately, 



S n ~ 0.44 + 0.22n (n > 



(80) 



We did not consider processes where a particle tunnels 
to a location other than a nearest-neighbor site, since the 
action for such processes will be considerably higher. The 
exchange frequency depends exponentially on the action, 
so that even a relatively small increase in S can suppress 
J quite substantially. 

We define the dimensionless prefactor A n by writing 




rV2 



s„ 



(81) 
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TABLE I. The dimensionless action S n for various ex- 
change processes, see Fig. 0. 



2 3 4 5 




6 6b 6c 6d 






SI) 




10 



11 



FIG. 2. The most important exchange paths (and some 
less important ones). The paths for n = 12, 14 and 16 can 
be found by adding a ring of particles around the n = 6, 8 
and 12 diagrams, in the same way as the n — 8, 9, 10 and 11 
diagrams can be derived from n — 2, 3, 4 and 6. 
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FIG. 3. The classical action for the first ten ring exchange 
processes. 
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FIG. 4. Scaling of the prefactor for two- to six-particle ring 
exchanges with inverse system size. N is the number of par- 
ticles that are allowed to move. 



where Ry = e 2 /2eaB is the effective Rydberg constant. 
In contrast to the classical action, the prefactor shows a 
strong dependence on the system size, as can be seen in 
Fig. 0. The TV-dependence fits well to a scaling form 



A„(N) 



Oi 

N 



a-1 

N 2 '' 



(82) 



from which we cam extract the values for the infinite sys- 
tem: A 2 = 2.60,B As = 2.19, i 4 = 2.48, A 5 = 3.15, and 
I 6 = 2.90. 



B. Results in the Presence of Disorder 



Although the technique outlined in Sec. Ill Bj allows 
us in principle to calculate the prefactor in the presence 
of impurities as well, we maintain the viewpoint that all 
qualitatively important changes in the exchange frequen- 
cies will be caused by variations of the classical action 
with disorder, and that the prefactor A n depends only 
weakly on disorder. This is by no means guaranteed, 
and in particular the dependence of the typical phonon 
frequencies on the disorder has to be investigated fur- 
ther. Note also that close to the melting transition an- 
harmonicities that are not accounted for in the instan- 
ton approximation may soften the phonon modes signif- 
icantly. We assume, however, that the disorder depen- 
dence of the exchange frequencies is dominated by the 
exponential factor. A small change in either the action 
or the prefactor causes a relative change 



SJn 



SA n 



^ /2 5, 



SS n 

On 



(83) 



in the exchange frequency. For moderately large values 
of r s the last term will give the dominant contribution. 
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FIG. 5. Two-particle exchange frequency relative to its 
value for the clean system, as a function of impurity layer 
distance 
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FIG. 7. Four-particle exchange frequency relative to its 
value for the clean system, as a function of impurity layer 
distance 
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FIG. 6. Three-particle exchange frequency relative to its 
value for the clean system, as a function of impurity layer 
distance 



To be explicit, we define the "reduced" exchange con- 
stants 



K„ 



(84) 



and study their dependence on disorder. Figs. |5j-|7] show 
the disorder averages of K n for n = 2, 3 and 4, nor- 
malized by their values K„ for the clean system. Here 
d/a is the distance to the impurity layer in units of the 
lattice constant of the clean Wigner crystal. The impu- 
rity concentration is taken to be x = 1/8 impurity ions 
per electron. The system size used for determining the 
equilibrium configuration is N = 48. 

While the impurities are practically of no effect for 
d/a > 1, in each case we see an enhancement of the 
average exchange frequency by up to a factor of 3 (at r s = 
30; at lower carrier densities the enhancement will be 



considerably larger) at d/a ~ 0.5. For smaller values of 
d/a three- and four-particle exchange frequencies decline 
again, and for d/a — > K-$ even falls below its original 
value. In the following we give an interpretation of these 
results; more details on the frequency distribution can be 
found in Appendix |Cj. 

In Fig. U we show the deviation of the electrons' equi- 
librium configuration from the ordered lattice, due to 
disorder. The average (static) root mean square dis- 
placement shows a sharp increase around the same value 
d/a ~ 0.5 at which the exchange frequencies peak. 
We claim that both signatures are due to a structural 
crossover that will be investigated in more detail in the 
following section. In the crossover region fluctuation ef- 
fects are amplified, which results in a softening of the 
potential barriers and an increase in the variance of the 
action of a instanton process. The increase in variance 
leads to an increase of the average exchange frequency, 
due to the positive curvature of J n (S). 

While three-and four-particle exchange frequencies fall 
off as we decrease d/a below 0.5, Ki remains enhanced 
by a factor of about 2.5 with respect to the clean system. 
Hence, if the distance to the impurity layer becomes very 
small, the two-particle exchange will dominate the mag- 
netic properties and enhance antiferromagnetic correla- 
tions. Of course antiferromagnetic order can be realized 
only locally. Nevertheless, this raises the possibility of 
a magnetic crossover signature that can, in principle, be 
picked up experimentally. 

Presumably this enhancement of the two-particle ex- 
change frequency is due to impurities mediating spin 
singlet correlations between the electrons of the Wigner 
glass (see also the next section). If an electron is trapped 
by an impurity charge, its repulsive interaction with 
neighboring electrons will be greatly reduced by the im- 
purity potential. Therefore exchange paths in which a 
second electron moves very close will not be suppressed 
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FIG. 8. Root mean square displacement of the electrons 
from their equilibrium positions, due to disorder. The three 
curves are for different system sizes. 



as much in the partition function. Unless the impurity 
concentration is very high (x ~ 1), this mechanism will 
not apply to n > 2 exchange processes. 



V. THE STRUCTURAL CROSSOVER 

For large d/a the influence of the impurities is weak, 
and on a local scale the lattice will be only slightly dis- 
torted (Fig. g). In the opposite limit (d/a — > 0), some 
electrons will be trapped in the potential wells created by 
impurity charges. These electrons are effectively removed 
from the lattice. The electron-impurity pairs now appear 
as dipoles with a dipole strength proportional to d, and 
therefore the effective disorder strength decreases as d 
becomes smaller. The remaining electrons will rearrange 
themselves to form a local Wigner lattice with an elec- 
tron density less than that of the clean system (Fig. 10). 
In the classical limit, and for d exactly equal to zero, the 
remaining electrons again form a perfect Wigner crystal, 
but with its electron density reduced by a factor I — x 1 
where x is the impurity concentration. 

The two limits d/a = oo and d/a — therefore corre- 
spond to two distinct structural phases of the system, 
with different translational symmetries. However, no 
phase transition involving these symmetries can occur 
at any finite value of d/a. It has been known for a long 
time that disorder destroys any spatial long-range order 
in two-dimensional systemsja and more recently it has 
been shown that quasi- long range— order, such as in a 
Bragg glass, cannot survive either .EJ To confirm the ab- 
sence of long range order we show in Fig. |llj the depen- 
dence of the rms deviations on the system size for weak 
disorder. The deviations increase linearly with the parti- 
cle number, hence quadratically in the linear dimensions. 

Locally, however, we can observe a sharp crossover 
between the two structural "phases". The correlation 
length will be strongly enhanced in the crossover region, 
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FIG. 9. Snapshot of the classical equilibrium electron con- 
figuration at d/a — 0.7. The filled circles show the electron 
positions, while the checked squares indicate the locations of 
the impurities. 
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FIG. 10. Snapshot of the classical equilibrium electron con- 
figuration at d/a = 0.1. The filled circles show the electron 
positions, while the checked squares indicate the locations of 
the impurities. 



but will remain finite due to a long distance cutoff im- 
posed by disorder. This enhancement of the correla- 
tion length will give rise to similar effects as are usually 
connected with phase transitions, such as the softening 
of phonon modes and enhancement of fluctuations. Of 
course these effects can only be observed locally. 



VI. CONCLUSIONS 

We have attempted to characterize the effective low en- 
ergy spin Hamiltonian for a disordered Wigner crystal, or 
a Wigner glass and have shown that disorder can make 
a qualitative difference. In particular, disorder can cause 
an enhancement of the two-particle exchange frequency 
relative to the other exchange frequencies, thereby mak- 
ing a possible ferromagnetically ordered state less likely 
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FIG. 11. rms deviation of electron coordinates from their 
equilibrium positions, due to disorder, at d/a — 2, for various 
system sizes. The dotted line is a least-square fit. 



and a spin liquid phase more likely. The solution of such 
a complex, competing multiparticlc spin Hamiltonian in- 
cluding disorder is a formidable many body problem. 
Nonetheless, it would be surprising if this Hamiltonian 
did not exhibit a multiplicity of competing phases in the 
ground state. We now present a speculative phase di- 
agram in the r s vs. disorder plane (see Fig. |l2| ), based 
primarily on symmetry considerations that can provide 
some guidance in the future. 

Let us first focus on the zero disorder axis: In the 
limit of very low densites (r s — > oo), three-particle ex- 
change will be most relevant, leading to a state with 
ferromagnetic long-range order (FM).H Upon increasing 
the density, two- and four-particle exchange will frus- 
trate the ferromagnet, until ferromagnetic order disap- 
pears at a critical density Tf. Within a mean-field ap- 
proximation, the (truncated) effective spin Hamiltonian 
can exhibit a variety _af multi-sublattice antiferromag- 
netic phases (MSAF).EZL-There is also the possibility of 
a spin liquid phase (SL).E3 For even higher densities, the 
Wigner crystal will quantum melt at r s = r c . 

While the ferromagnet with disorder averaged order 
parameter < S, > 7^ can survive in the presence of 
weak disorder, the various MSAF states are tied to the 
existence of translational symmetry. As noted earlier, 
crystalline order will be destroyed on long length scales 
by any amount of disorder in two dimensions. Hence, 
no multi-sublattice antiferromagnetic phases can exist in 
the disordered system. Similarly, various dimerized bro- 
ken symmetry states cannot be distinct states of matter 
either .l 2 ] We expect a crossover region, with short-range 
antiferromagnetic correlations, to a spin glass phase (SG) 
at T = 0, although there may not be a finite temperature 
spin glass transition in d = 2. (In fact, numerical finite- 
size diagonalizations in simpler randomly frustrated spin- 
1/2 nearest neighbor Heisenberg models have been ai- 
gued to exhibit spin glass behavior in the ground state.cij) 
The complexity of this problem can be visualized by the 




•s 

FIG. 12. Speculative phase diagram in the r s - disorder 
plane. The various phases are: Anderson insulator (AI), para- 
magnetic (Efros-Shklovskii) insulator (PI), spin glass (SG), 
spin liquid (SL), various multi-sublattice antiferromagnetic 
states (MSAF), ferromagnet (FM), and a metallic state with 
broken time-reversal symmetry (M). The diagonally shaded 
area is a crossover region with intermediate-range MSAF 
correlations. S and T indicate broken spin-rotational and 
time-reversal symmetry, respectively. Not included is a pos- 
sible superconducting state with broken U(l) symmetry. 



fact that quenched disorder in such a frustrated multipar- 
ticle exchange Hamiltonain at T = appears as infinitely 
long ranged correlated disorder in the imaginary time di- 
rection in the field theoretical description. In any case, 
the spin glass phase should be characterized by a nonva- 
nishing Edwards-Anderson order parameter, that is, the 
disorder average < S t > 2 is nonzero, but < Sj > = 0. 

It is also interesting to construct a chiral order param- 
eter. Let R4, R 2 , and R3 be the vertices of a triangular 
plaquette and define $(R) = (S(Ri) • S(R 2 ) x S(R 3 )), 
where R is the center of the triangle. In the spin glass 
phase, <&(R) = 0, but parasitically $ 2 (R) 7^ 0. This 
suggests a transition to an adjacent phase in which the 
Edwards- Anderson order parameter is zero, but $ 2 (R) 7^ 
0. This is a new state of matter with broken time rever- 
sal symmetry, which we shall label to be the random flux 
state. It is interesting to ask what a prototype model 
could be. It is tempting to speculate that this is a vari- 
ant of the random flux model. Although there is some 
evidence of a metal-insulator transition in the random 
flux model, separating a T-broken metallic state from a 
T-broken insulating state, the evidence to the contrary 
also existso The resolution of this controversy should be 
an important advance. 

Next, we look at the strong-disorder, low-density re- 
gion. Since the exchange frequencies decrease exponen- 

1/2 
tially with r s , the characteristic energy scale of disor- 
der will be the largest energy scale, so that the carriers 
are independently trapped in some minima of the dis- 
order potential. The resulting state will be a paramag- 
netic (Efros-Shklovskii) insulator. Since no symmetries 
are broken in this state, it can be continuously connected 
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M 

T~ 

4 
6 
8 

12 
16 
32 



S 3 
1.1924 
1.4919 
1.5455 
1.5638 
1.5763 
1.5804 
1.5840 



-25% 
-5.9% 
-2.5% 
-1.3% 
-.54% 
-.28% 
-.05% 



TABLE II. Dependence of the classical action on the num- 
ber of time slices used in the calculation. 



to the noninteracting disordered electron system, which 
is an Anderson insulator. ._. 

A superconducting state with broken U(l) symmetryEEl 
is possible in principle. It is not clear to us, however, 
where in the phase diagram such a superconducting state 
should occur, therefore it is not included in the phase 
diagram shown in Fig. (12). We certainly do not imply 
that such a phase is impossible. 

In none of the phases involving broken T and 5, the 
impurity potential can couple to the order parameter as a 
"random field" . Rather, the effect of the potential scat- 
tering due to impurities is to randomize the exchange 
constants. Similar to_dgorous results known in classical 
statistical mechanics,o we can argue that these broken 
symmetry transitions in the ground state are necessarily 
continuousHj Thus, scaling must hold at these quantum 
phase transitions, and the signature of quantum critical- 
ity should be observable at finite, but low temperatures. 
In contrast, the Wigner crystal transition of a pure sys- 
tem is a first order transition at which scaling could not 
possibly hold. 
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APPENDIX A: ACCURACY ISSUES 

In Table || we see how the classical action depends on 
the number of time slices M used in the discretized form 
(p7|). We use the three-particle exchange, with 27 mobile 
particles, as an example. The errors for other exchange 
processes are comparable. Extrapolating to M = oo 
yields S3 ~ 1.5848. Errors are given with respect to 
this value. 



-^mobile 


L 


S3 


err 


3 





1.93334 


24% 


12 


1 


1.65910 


6.6% 


27 


2 


1.57627 


1.3% 


46 


3 


1.56166 


.38% 


69 


4 


1.55762 


.11% 


96 


5 


1.55628 


.03% 



TABLE III. Dependence of the classical action on the num- 
ber of electrons allowed to move. 



Table III shows the dependence on the number of par- 
ticles allowed to move in the exchange. Here L is the 
number of layers added around the interchanging parti- 
cles. The number of time slices is M = 12. Extrapolation 
to L = 00 yields S3 ~ 1.5558. Since the errors introduced 
by finite M and finite ./Vmobiic are of opposite sign, the 
choices M = 16 and N moh \x c ~ 80 for the clean system 
should yield accurate results to within ~ 0.2 percent. 
For the disordered system, where we used M = 8 and 
-^mobile — 32, we expect the systematic errors to be on 
the order of 1 percent. 



APPENDIX B: PREFACTOR IN ONE 
DIMENSION 



Here we apply the method introduced in Sec. Ill B 
for numerical calculation of the prefactor to tunneling 
in a double well potential in one dimension, which can 
be solved exactly. This serves as a useful test for the va- 
lidity of the technique. The double well potential is given 
by 



V(x) = {x 2 - l) 2 



(Bl) 



and the mass of the particle is set to m = 1 . According 
to Ref. pi a general expression for the prefactor P in one 
dimension is 



P- 



-uj T t 



A ' 



(B2) 



where T T is the length of the time slice (we set T T — > 00 
at the end of the calculation) and ojq = \Z(PV/dx 2 \ x —i 
is the attempt frequency. Ao is the lowest eigenvalue of 
the equation 



[-d* + V(x(T))}u(T)=Au(T), 

which is given byc3 

ALUn 



An 



"°_ A 3 e -a*T 



(B3) 



(B4) 



where S; nst = J dx y f 2V{x) is the action along the classi- 
cal trajectory, and A is the prefactor for the asymptotic 
form of the first time derivative of the classical trajectory: 
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M 


4 


8 


16 


32 


64 


128 


A 


10.25 


9.970 


9.859 


9.819 


9.804 


9.799 



TABLE IV. Prefactor in the double well problem for dif- 
ferent numbers of time slices. 



d/a 


s 2 


02 


s 3 


0-3 


Si 


(74 


2.0 


1.631 


0.025 


1.521 


0.016 


1.662 


0.023 


1.5 


1.633 


0.043 


1.519 


0.031 


1.659 


0.042 


1.0 


1.627 


0.089 


1.514 


0.066 


1.649 


0.083 


0.7 


1.621 


0.144 


1.501 


0.120 


1.626 


0.138 


0.6 


1.614 


0.171 


1.490 


0.153 


1.600 


0.177 


0.5 


1.587 


0.210 


1.487 


0.185 


1.588 


0.231 


0.4 


1.586 


0.227 


1.513 


0.184 


1.627 


0.258 


0.3 


1.572 


0.234 


1.533 


0.189 


1.678 


0.290 


0.2 


1.598 


0.254 


1.569 


0.189 


1.725 


0.297 



TABLE V. Dependence of mean and standard deviation of 
the dimensionless action on disorder strength. 

± c (t) ~ Ae~" oM as t -> ±oo. (B5) 

Integrating the equation of motion ( |l6| ) yields 

x c (t) = tanhV2r. (B6) 



From Eqs. (B2), (B4) and (B6) we immediately get 

P = 4^6 ~ 9.798. (B7) 



Table IV shows results of a numerical computation us- 



ing M time slices. We see that the technique reproduces 
the exact result within one percent accuracy for M > 16. 



APPENDIX C: DISTRIBUTION OF EXCHANGE 
FREQUENCIES 

The dimensionless action S n that enters the expression 
( p4| ) for the exchange frequencies depends on the particu- 
lar realization of disorder, and can therefore be viewed as 
a random variable. In most regions of parameter space 
the random distribution turns out to be well-described 
by a normal distribution 



P(Sn) 



1 



(gn- 



2ir<7„ 



(CI) 



In this approximation, the frequency distribution of 
S n , and thereby of K n , can be reconstructed from two 
parameters, the mean action S n and the standard devia- 
tion cr„, which are listed in Table |v| for various disorder 
strengths. 

For d/a ^0.1 the random distribution acquires a sig- 
nificant non-gaussian component, hence the correspond- 
ing values are not listed. 
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